{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Supervised Learning with Neural Networks\n",
    "\n",
    "In this tutorial we show how to optimize a neural networks to approximate the state given. In this example, we consider the ground state of the J1-J2 model in one-dimension obtained from ED using NetKet. \n",
    "\n",
    "The Hamiltonian of the model is given by:\n",
    "$$ H = \\sum_{i=1}^{L} J_{1}\\hat{S}_{i} \\cdot \\hat{S}_{i+1} + J_{2} \\hat{S}_{i} \\cdot \\hat{S}_{i+2} $$\n",
    "where the sum is over sites of the 1-D chain.\n",
    "\n",
    "\n",
    "## Outline:\n",
    "    1. Obtain data from ED\n",
    "    2. Choosing the machine (variational ansatz) and the optimizer\n",
    "    3. Defining the Supervised Learning object\n",
    "    4. Running the Supervised Learning\n",
    "    5. Data Visualisation\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import netket library\n",
    "import netket as nk\n",
    "\n",
    "# Helper libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Obtain data from ED\n",
    "\n",
    "For a supervised learning problem, we would need to provide the data including input $X$ and output label $Y$. The neural network is asked to learn the mapping $Y=f(X)$. In our case, the input is the spin basis and the output is the coefficient of the corresponding spin basis. \n",
    "\n",
    "First, we write a simple function to obtain data, i.e. ground state, from exact diagonalization. For detailed explanation see the tutorial for J1-J2 model for example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_ed_data(L, J2=0.4):\n",
    "    # Sigma^z*Sigma^z interactions\n",
    "    sigmaz = np.array([[1, 0], [0, -1]])\n",
    "    mszsz = (np.kron(sigmaz, sigmaz))\n",
    "\n",
    "    # Exchange interactions\n",
    "    exchange = np.asarray(\n",
    "        [[0, 0, 0, 0], [0, 0, 2, 0], [0, 2, 0, 0], [0, 0, 0, 0]])\n",
    "\n",
    "    # Couplings J1 and J2\n",
    "    J = [1., J2]\n",
    "\n",
    "    mats = []\n",
    "    sites = []\n",
    "\n",
    "    for i in range(L):\n",
    "        for d in [0, 1]:\n",
    "            # \\sum_i J*sigma^z(i)*sigma^z(i+d)\n",
    "            mats.append((J[d] * mszsz).tolist())\n",
    "            sites.append([i, (i + d + 1) % L])\n",
    "\n",
    "            # \\sum_i J*(sigma^x(i)*sigma^x(i+d) + sigma^y(i)*sigma^y(i+d))\n",
    "            mats.append(((-1.)**(d + 1) * J[d] * exchange).tolist())\n",
    "            sites.append([i, (i + d + 1) % L])\n",
    "\n",
    "    # 1D Lattice\n",
    "    g = nk.graph.Hypercube(length=L, n_dim=1, pbc=True)\n",
    "\n",
    "    # Spin based Hilbert Space\n",
    "    hi = nk.hilbert.Spin(s=0.5, graph=g)\n",
    "\n",
    "    # Custom Hamiltonian operator\n",
    "    ha = nk.operator.LocalOperator(hi)\n",
    "    for mat, site in zip(mats, sites):\n",
    "        ha += nk.operator.LocalOperator(hi, mat, site)\n",
    "\n",
    "    # Perform Lanczos Exact Diagonalization to get lowest three eigenvalues\n",
    "    res = nk.exact.lanczos_ed(ha, first_n=3, compute_eigenvectors=True)\n",
    "\n",
    "    # Eigenvector\n",
    "    ttargets = []\n",
    "\n",
    "    tsamples = []\n",
    "\n",
    "    for i, visible in enumerate(hi.states()):\n",
    "        # only pick zero-magnetization states\n",
    "        mag = np.sum(visible)\n",
    "        if(np.abs(mag) < 1.0e-4):\n",
    "            tsamples.append(visible.tolist())\n",
    "            ttargets.append([np.log(res.eigenvectors[0][i])])\n",
    "\n",
    "    return hi, tsamples, ttargets\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After we obtain the result as ``res``, we return the hilbert space ``hi``, the spin basis ``tsamples``, and the coefficients ``ttargets``.\n",
    "\n",
    "Notice that we restrict ourselves to $\\sum S_z = 0$ symmetry sector to simplify the learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now consider a small system $L=10$ and with $J_2 = 0.4$, and obtain the data by calling the function ```load_ed_data```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 10\n",
    "J2 = 0.4\n",
    "\n",
    "# Load the Hilbert space info and data\n",
    "hi, training_samples, training_targets = load_ed_data(L, J2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) Choosing the Machine and the Optimizer\n",
    "\n",
    "For this tutorial, we consider the Restricted Bolzmann Machine ``nk.machine.RbmSpin`` and the AdaDelta optimizer ``nk.optimizer.AdaDelta``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Machine\n",
    "ma = nk.machine.RbmSpin(hilbert=hi, alpha=1)\n",
    "ma.init_random_parameters(seed=1234, sigma=0.01)\n",
    "# Optimizer\n",
    "op = nk.optimizer.AdaDelta()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Defining the Supervised Learning object\n",
    "\n",
    "We have now have almost everything (machine, optimizer, data) for setting up a supervised learning object. We also need to provide the batch size, ``batch_size``, for stochatic gradient descent. For detail, see https://en.wikipedia.org/wiki/Stochastic_gradient_descent\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Supervised learning object\n",
    "spvsd = nk.supervised.Supervised(\n",
    "    machine=ma,\n",
    "    optimizer=op,\n",
    "    batch_size=400,\n",
    "    samples=training_samples,\n",
    "    targets=training_targets)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4) Running the Supervised Learning\n",
    "\n",
    "The very last piece we need for supervised learning is the loss function.\n",
    "\n",
    "### Loss function\n",
    "There are different loss functions one could define for the optimization problem, for example:\n",
    "$$\n",
    "\\begin{align*}\n",
    "        \\mathcal{L}_\\text{MSE log} &= \\frac{1}{N} \\sum_{i}^N |\\log\\Psi(X_i) - \\log\\Phi(X_i) |^2\\\\\n",
    "        \\mathcal{L}_\\text{Overlap} &=-\\log\\Big[ \\frac{\\langle{\\Psi|\\Phi}\\rangle\\langle{\\Phi|\\Psi}\\rangle}{\\langle{\\Psi|\\Psi}\\rangle\\langle{\\Phi|\\Phi}\\rangle} \\Big] \\\\\n",
    "        &=- \\log\\Big( \\sum_{i}^N \\Psi^*(X_i)\\Phi(X_i) \\Big) - \\log\\Big( \\sum_{i}^N \\Phi^*(X_i)\\Psi(X_i) \\Big) \\\\\n",
    "        &\\qquad +\n",
    "        \\log\\Big( \\sum_{i}^N \\Psi^*(X_i)\\Psi(X_i) \\Big) +\n",
    "        \\log\\Big( \\sum_{i}^N \\Phi^*(X_i)\\Phi(X_i) \\Big)\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Here, we consider the latter one, which is the negative log of the overlap, as the loss function.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "### Gradient estimate\n",
    "Taking the derivative from overlap errror function above, we have\n",
    "$$\n",
    "\\begin{equation*}\n",
    "    \\partial_k \\mathcal{L}_\\text{Overlap} = -\\frac{\\sum_i O_k^*\\Psi^*(X_i)\\Phi(X_i) }{\\sum_i\\Psi^*(X_i)\\Phi(X_i)} + \\frac{\\sum_i O_k^*\\Psi^*(X_i)\\Psi(X_i)}{\\sum_i \\Psi^*(X_i)\\Psi(X_i)}\n",
    "\\end{equation*}\n",
    "$$\n",
    "\n",
    "Note that $N$ is the size of the Hilbert space. In general, this could not be computed exactly.\n",
    "\n",
    "\n",
    "We could estimate this gradient by sampling different distributions,\n",
    "$$\n",
    "\\begin{equation*}\n",
    "    \\hat{\\partial_k \\mathcal{L}}_\\text{Overlap uni} = \\frac{\\Big\\langle O_k^*\\Psi^*(X_i)\\Psi(X_i)\\Big \\rangle_{i\\sim\\text{uni}[1,N]} }{\\Big \\langle \\Psi^*(X_i)\\Psi(X_i) \\Big \\rangle_{i\\sim\\text{uni}[1,N]}} - \\frac{\\Big \\langle O_k^*\\Psi^*(X_i)\\Phi(X_i)\\Big \\rangle_{i\\sim\\text{uni}[1,N]} }{\\Big \\langle \\Psi^*(X_i)\\Phi(X_i) \\Big \\rangle_{i\\sim\\text{uni}[1,N]}} \n",
    "\\end{equation*}\n",
    "$$\n",
    "$$\n",
    "\\begin{equation*}\n",
    "    \\hat{\\partial_k \\mathcal{L}}_\\text{Overlap phi} = \\frac{\\Big \\langle O_k^*(X_i)\\frac{\\lVert \\Psi(X_i)\\rVert^2}{\\lVert \\Phi(X_i)\\rVert^2} \\Big \\rangle_{X_i\\sim \\lVert \\Phi(X_i)\\rVert^2 }}   {\\Big \\langle \\frac{\\lVert \\Psi(X_i)\\rVert^2}{\\lVert \\Phi(X_i)\\rVert^2} \\Big \\rangle_{X_i\\sim \\lVert \\Phi(X_i)\\rVert^2 }} - \\frac{\\Big \\langle O_k^*(X_i)\\frac{ \\Psi^*(X_i)}{ \\Phi^*(X_i)} \\Big \\rangle_{X_i\\sim \\lVert \\Phi(X_i)\\rVert^2 }}{\\Big \\langle \\frac{ \\Psi^*(X_i)}{ \\Phi^*(X_i)} \\Big \\rangle_{X_i\\sim \\lVert \\Phi(X_i)\\rVert^2 }}\n",
    "\\end{equation*} \n",
    "$$\n",
    "\n",
    "\n",
    "\n",
    "So for the overlap loss function, we have two gradient estimate, one is $\\hat{\\partial_k \\mathcal{L}}_\\text{Overlap uni}$, ```Overlap_uni```, and $\\hat{\\partial_k \\mathcal{L}}_\\text{Overlap phi}$, ```Overlap_phi```."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We save the loss function every iteration, and save the optimized parameters only every ``save_params_every`` iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of iteration\n",
    "n_iter = 4000\n",
    "\n",
    "# Run with \"Overlap_phi\" loss. Also available currently is \"MSE, Overlap_uni\"\n",
    "spvsd.run(n_iter=n_iter, loss_function=\"Overlap_phi\",\n",
    "          output_prefix='output', save_params_every=50)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5) Data Visualisation\n",
    "\n",
    "We have optimized our machine to approximate the ground state of the J1-J2 model. The results for the loss function are stored in the \".log\" file and the optimized parameters in the \".wf\" file. The files are all in json format.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data from the .log file\n",
    "import json\n",
    "\n",
    "data=json.load(open(\"output.log\"))\n",
    "\n",
    "# Extract the relevant information\n",
    "iters=[]\n",
    "log_overlap=[]\n",
    "mse=[]\n",
    "mse_log=[]\n",
    "\n",
    "data=json.load(open('output.log'))\n",
    "for iteration in data[\"Output\"]:\n",
    "    iters.append(iteration[\"Iteration\"])\n",
    "    log_overlap.append(iteration[\"log_overlap\"])\n",
    "    mse.append(iteration[\"mse\"])\n",
    "    mse_log.append(iteration[\"mse_log\"])\n",
    "\n",
    "overlap = np.exp(-np.array(log_overlap))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we plot the overlap, i.e. fidelity, with respect to the number of iteration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEYCAYAAABLOxEiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl4U1X6wPHv25ZS1kLZ94KAyL6LIioCCoIi6rg77uig4za/cXDQGbeZUUdFcVBHR8V9ARVRUEQFQVFW2QVa9rKVtYUW6Pb+/shtSUvaJm2Sm5b38zx5cnNyc++bmzZv7jnnniOqijHGGOOvKLcDMMYYU7FY4jDGGBMQSxzGGGMCYonDGGNMQCxxGGOMCYglDmOMMQGxxGGMMSYgljiMMcYExBKHMUEkIpNE5Ak/190sIoNDHZMxwWaJw0QkEWkoIioijYOxXkUX6vcpIgki8pmIZIjIFhG5xo/XtBORoyLybihiMpErxu0AjClGN2CPqu4K0noVXajf50QgC2gEdAemi8hyVV1dymsWhSgeE8HsjMNEqq7AiiCuV8CpIvqziKxwfmG/LiKNROQrETkkIt+KSF1n3dNEZI6IHBSR1SJycZFt9RCRpc7rPgLivJ5rKiKfiMgeEdkkIncHEmd536e/RKQGcBnwsKoeVtUfgWnA9SW85irgIPBdKGIykc0Sh4lUXfDvi9Lf9Yq6DBgCtAcuAr4C/grUx/N/cbeIVAG+AL4BGgJ/BN4TkVMBRCQWmAq8AyQAk53tIiJRzmuXA82AQcC9InJBGWKFAN6niHzpJDpfty99vKQ9kKuq673KlgOditl+beAx4E8BvgdTSVjiMJGq4Be2iMSLyEIROSwinX2tJyJniMjPIvKDiHzgfOmX5EVV3a2q24F5wAJV/VVVjwGfAT2AfkBN4ElVzVLV74EvgaudbfQDqgDPq2q2qk7heNVNH6CBqj7mvHYj8BpwVXmOhz/vU1VHqGqdYm4jfGy7JpBWpCwNqFVMLI8Dr6vqtjK+F1PBWeIwEUdEooHT8PzqBcgEhgNTSlhvC3Ceqp4DbARGlrKb3V7LR3w8rgk0Bbapap7Xc1vwnEHgPL9dC89NsMW5bwU09f61j+eMplEpcZ2gnO/TH4eB2kXKagOHfMTSHRgMjA/Cfk0FZY3jJhK1x/O3uQZAVbOBPSJS7HrOmUK+HCCv6MplsANoISJRXsmjJZBfpbMTaCYi4pU8WgIbgG3AJlVtF4Q4AnqfIvIVMKCYbc1T1WFFytYDMSLSTlWTnLJugK+G8XOBRGCr83nUBKJFpKOq9vTv7ZiKzs44TCTqCqwv8iXp13oi0hoYhqdKqbwWABnAAyJSRUTOxdMe8qHz/M94vrzvFpEYEbkU6Os8txBIF5G/iEg1EYkWkc4i0sfXjpzrPyYVE0dA71NVh6lqzWJuRZMGqpoBfAo8JiI1RKQ/njOZd3zE8ipwCp6eV92BV4DpQFnbbkwFZInDRKIuHK+m8ns9p9H2LeB6Vc0qbxDONi7G8wW9F3gJ+L2qrvV6/lLgRuAAcCWeL2BUNRdPkukObHJe/z8gvpjdtQB+Kua5kL5PxxigGpAKfAD8Ib8rrtPb7K8Aqpqpqrvyb3iquY6q6p4gxWEqALGpY02kEZHZwPuq+lqR8knAM6q6quh6IhIDfA486zRiVxhO76zlQFenWq7o85XifZrKwxKHiSgiMgTPL94OqrrXq3wGnl/vW4D/Atu91xOR6/E02K5yXvKyqn4U1uBDoOjxqKzv01QsljhMxBCRlUAucK+qzinvehXdyfI+TcVjicMYY0xArHHcGGNMQCrldRz169fXxMREt8MwxpgKZcmSJXtVtUFp61XKxJGYmMjixYvdDsMYYyoUEdlS+louV1WJyBsikioiq4p5XkRkgogkOyOZ2pWpxhjjMrfbOCYBQ0t4fhjQzrmNBl4OQ0zGGGNK4GpVlarOFZHEElYZCbztjAP0i4jUEZEmqrqzpO0uWbIEH+MaGWOMCQK3zzhK0wzPYHH5Ujg+MmkhIjJaRBaLiDVuGGNMCEV647iv0wafF56o6qt4BmCjd+/eao3jkSPlQCZLtx5kxbaDLE85yKrt6RzJzi14Pq5KFA1rxdGwVlUa1q5KneqxVK8STfWqMdSI9dxXrxJNjarRVI+NoUbVaKpViSn0OC4mmqgoO8s0pjz8ramJ9MSRgmfwt3zN8Qx1bSLY0excftm4jx/W7+GH9XvYuCcDgNiYKDo3rc2VfVrQpVk8bRvWJLFeDWpXi7GqRWMqkEhPHNOAu0TkQ+B0IK209g3jDlVl4ab9TF6SwoyVO8nMyqVqTBSnt6nHNX1b0q9NPU5tXIsq0ZFeO2qMKY2riUNEPsAzMUx9EUkB/o5nKk5U9RVgBnAhkIxnFrib3InUFGfHwSN8siSFKUtT2LIvkxqx0Yzo2oRhXZpwRpt6xFWJdjtEY0yQud2r6upSnlfgzjCFYwKwZV8GL36fzKdLU8hTOKNNPe4+rx3DujSmemykn8gaU3nNnTuXe++9lxUrVvDhhx9y+eWXB30f9h9uAvLbznRemrOBGSt3EhMl3NS/NTeemUiLhOpuh2aM63Jzc4mOdvcsu2XLlkyaNIlnnnkmZPuwCmfjl9T0o/xlygounDCPOWtTueWs1sx7YCAPj+hoScOU2ebNm+nQoQO33nornTt35tprr+Xbb7+lf//+tGvXjoULFwKwcOFCzjzzTHr06MGZZ57JunXrAHjuuee4+eabAVi5ciWdO3cmMzPzhH0MGDCAnj170rNnT+bPn1/w3NNPP02XLl3o1q0bY8eOBSA5OZnBgwfTrVs3evbsyYYNG5gzZw4jRowoeN1dd93FpEmTAM8QR4899hhnnXUWkydP5rXXXqNPnz5069aNyy67rCCe3bt3M2rUKLp160a3bt2YP38+Dz/8MC+88ELBdseNG8eECRPKdUwTExPp2rUrUVEh/HpX1Up369Wrl5rgyDiWreNnrdPTHv5K2/51uj7+xWo9mJHldlgmBPB0dQ/6rSSbNm3S6OhoXbFihebm5mrPnj31pptu0ry8PJ06daqOHDlSVVXT0tI0OztbVVVnzZqll156qaqq5ubm6oABA/TTTz/VXr166Y8//njCPjIyMvTIkSOqqrp+/XrN/36YMWOGnnHGGZqRkaGqqvv27VNV1b59++qnn36qqqpHjhzRjIwMnT17tg4fPrxgm3feeae++eabqqraqlUrfeqppwqe27t3b8HyuHHjdMKECaqqesUVV+j48eNVVTUnJ0cPHjyomzZt0h49ehS8lzZt2hR6fb6zzjpLu3XrdsJt1qxZxR7bG264QSdPnlzs874Ai9WP71irqjLFWrsrnTveWcLmfZkM69yYB4Z2oHX9Gm6HZSqZ1q1b06VLFwA6derEoEGDEBG6dOnC5s2bAUhLS+OGG24gKSkJESE72zPDblRUFJMmTaJr167cfvvt9O/f/4TtZ2dnc9ddd7Fs2TKio6NZv349AN9++y033XQT1at7zpgTEhI4dOgQ27dvZ9SoUQDExcX59R6uvPLKguVVq1bx0EMPcfDgQQ4fPswFF1wAwPfff8/bb78NQHR0NPHx8cTHx1OvXj1+/fVXdu/eTY8ePahXr94J2583b55fcYSLJQ5zguzcPF6Zs4EJ3ycRX60KH9zWjzNOOfGP2VQu6tKkblWrVi1YjoqKKngcFRVFTk4OAA8//DADBw7ks88+Y/PmzZx77rkFr0lKSqJmzZrs2OH7Eq/x48fTqFEjli9fTl5eXkEyUNUTrh8q7hjExMSQl5dX8Pjo0aOFnq9R4/gPqhtvvJGpU6fSrVs3Jk2axJw5c0p8/7feeiuTJk1i165dBdVuRQ0YMIBDhw6dUP7MM88wePDgErcfCtbGYQr5bWc6l0z8iWdnrWdo5yZ8c985ljSM69LS0mjWzDPaUH7bQn75Pffcw9y5c9m3bx9Tpkzx+domTZoQFRXFO++8Q26uZ9SC888/nzfeeKOgDWL//v3Url2b5s2bM3XqVACOHTtGZmYmrVq1Ys2aNRw7doy0tDS+++67YmM9dOgQTZo0ITs7m/fee6+gfNCgQbz8smec1tzcXNLT0wEYNWoUX3/9NYsWLSo4Oylq3rx5LFu27ISbG0kDLHEYx5GsXJ7+ei0Xvfgju9OP8sp1PXnx6h4k1Ih1OzRjeOCBB3jwwQfp379/wRc/wH333ceYMWNo3749r7/+OmPHjiU1NbXQa8eMGcNbb71Fv379WL9+fcHZwdChQ7n44ovp3bs33bt3L+iF9M477zBhwgS6du3KmWeeya5du2jRogVXXHEFXbt25dprr6VHjx7Fxvr4449z+umnM2TIEDp06FBQ/sILLzB79my6dOlCr169WL16NQCxsbEMHDiQK664Iig9shYtWkTz5s2ZPHkyt99+O506dSr3NouqlHOO21hVgZmzLpWHP1/Ftv1HuLxXc8ZdeBp1LWEYExZ5eXn07NmTyZMn065dO1djEZElqtq7tPWsjeMktj8ji0e/WM3ny3bQpkENa8swJszWrFnDiBEjGDVqlOtJIxCWOE5CuXnKZ79u558zfiP9SDb3DGrHmIGnUDXGhgcxJpw6duzIxo0b3Q4jYJY4TiIZx3KYtnwHr/+4ieTUw3RrUYenLutCh8a13Q7NGFOBlCtxiMg3qnq+s/ygqv4rOGGZYNqdfpS3f97MW/O3cPhYDh0a1+Kla3sytFNjm8PCGBOw8p5xNPBa/h1giSOCrN99iP/+sJHPl20nJ08Z3rUJN/dvTc+WdWz+C2NMmZU3cVS+LlkVnKqyaPMB/vvDBr5bm0pclSiu69eKG89MJNGu+jbGBEF5E0cbEZmGZ4rX/OUCqnpxObdv/LTn0DG+WbOLT5aksHTrQepWr8K9g9vx+zMS7VoMY0xQlTdxjPRaDt0YvsanI1m5TF+5k0+WpLBw835y85RW9arz6MWduKJ3C6rFWi8pY0zwlStxqOoPwQrE+EdVWbbtINOW7+CzX7dzMDObNvVrcMc5bbioW1NObVTL2i+MMSHlendcERkKvABEA/9T1SeLPN8SeAuo46wzVlVnhD1QlyWnHmbqr9uZtnwHW/dnEhsTxeDTGnLd6a3o16ae9Y4yxoSN23OORwMTgSFACrBIRKap6hqv1R4CPlbVl0WkI555yBPDHmyYqSpb9mXyy8Z9TF+5k3lJe4kS6N+2Pn88ry0XdG5M7bgqbodpjDkJuX3G0RdIVtWNACLyIZ52E+/EoUD+FWrxgO+xkyuBPYeOsWjzfpY7VVE70zxDNzerU437h7Tn6r4taVCrailbMcaY0Apa4hCRB1T16fx7P1/WDNjm9TgFOL3IOo8A34jIH4EagM9xhEVkNDAaPHPuVhSHj+Xw4ndJLN5ygCVbDgAQGx1F9xZ1+ON57ejZqo61WxhjIkowzziuAp72uveHr2/DoteGXA1MUtVnReQM4B0R6ayqeYVepPoq8Cp4RscNKHIX5OYpb83fzKT5m9m6P5NuzeO5f0h7BrSrT8emtW3cKGNMxApFVVUgP41TgBZej5tzYlXULcBQAFX9WUTigPpAKhWQqjJn3R7GfbaSHWlH6dC4Fi9f25NhXZq4HZoxxvjF7TaORUA7EWkNbMdztnJNkXW2AoOASSJyGhAH7AlrlEEyP3kvE+ck81PyPurXjOWJSzpz7ektrRrKGFOhuJo4VDVHRO4CZuLpavuGqq4WkceAxao6DfgT8JqI3IenGutGrUCzTx3JyuW1eRv5YvkOklIPE1clirvPa8ud57W16ihjTIXk9hkHzjUZM4qU/c1reQ3QP9xxBcPR7Fx+99/5rNqeTt/EBO4a2JbLejWntY0ZZYypwIKZOOY497ODuM0K7f6Pl7Fqezrjr+zGqB7N3Q7HGGOCImiJQ1Xv974/meXlKX+btooZK3fxpyHtLWkYYyoV16uqKhtV5a4PljJj5S4u6taUMQPbuh2SMcYElSWOIJu+ciczVu5iQLv6TLiqu/WYMsZUOlFuB1CZpGVm8/fPV5NYrzpv3tjHkoYxplIK5pAjlwJn4eky+6OqfhasbVcUk5dsY19GFpNu6ktMtOVkY0zlFJRvNxF5CbgDWAmsAm4XkYnB2HZFsXHPYZ6btZ4+iXXp0jze7XCMMSZkgnXGcQ7QOf/CPBF5C08SOWmM/XQlAoy/srvboRhjTEgFqz5lHeA9JG0LYEWQth3xklMPs3DTfm45qzXN61Z3OxxjjAmpYJ1x1AN+E5GFzuM+wM8iMg1AVS8O0n4i0sTZyQBc1bfiDOdujDFlFazE8bfSV6m8Fm7az4VdGtO0TjW3QzHGmJALSuJQ1R+CsZ2KKOVAJtsPHuG2Aa3dDsUYY8IiWL2q+onIIhE5LCJZIpIrIunB2Hakm7IkBYAz29Z3ORJjjAmPYDWO/wfPTH1JQDXgVqesUjuYmcXz3yYxoF192jeq5XY4xhgTFsEc5DBZRKJVNRd4U0TmB2vbkeq5WesBuHVAG5cjMcaY8AlW4sgUkVhgmYg8DewEKvWkE6rKjJU7GdqpMee0b+B2OMYYEzbBqqq6Hs8MfncBGXiu47gsSNuOSOt2H2Lv4SzOamdtG8aYk0uwelVtcRaPAI8GY5uR7oVvkwC4oFNjlyMxxpjwKlfiEJGVeAY19ElVu/qxjaHAC3jOWP6nqk/6WOcK4BFnX8tV9ZqyxhwMBzKy+H5tKt1a1KFBrapuhmKMMWFX3jOOEeV5sYhEAxOBIUAKsEhEpjnzjOev0w54EOivqgdEpGF59hkMSamHOZaTxz2DbJImY8zJp1yJQ1W3OF/+M1V1cBk20RdIVtWNACLyITASWOO1zm3ARFU94OwztTwxB8OPyXuJEujczEbBNcacfMrdOO50v80UkbJ8izYDtnk9TnHKvLUH2ovITyLyi1O1dQIRGS0ii0Vk8Z49e8oQin+ycvKY8F0SXZvXoWGtuJDtxxhjIlWwuuMeBVaKyCw8vaoAUNW7S3mdrynyiraZxADtgHOB5sA8EemsqgcLvUj1VeBVgN69exfb7lJeHy7aCsDg01yvMTPGGFcEK3FMd26BSsHTdTdfc2CHj3V+UdVsYJOIrMOTSBaVJdDy+jFpL63qVeeu89q5sXtjjHFdsLrjviUi1YCWqrougJcuAtqJSGtgO3AVULTH1FQ8w5lMEpH6eKquNgYh7IAdOprNgk37Oa+DnW0YY05ewRrk8CJgGfC187h7/lwcJVHVHDwXDc4EfgM+VtXVIvKYiOTP4TET2Ccia4DZwJ9VdV8w4g7UE1/+RtqRbG48M9GN3RtjTEQIVlXVI3h6SM0BUNVlzllEqVR1BjCjSNnfvJYVuN+5uWb97kN8tHgbdapXoVuLOm6GYowxrgrWkCM5qppWpCxkDdRu+GChp1H8jRv7uByJMca4K1hnHKtE5Bog2rlg726g0o2OGxsdRc+Wdd0OwxhjXBWsM44/Ap2AY8D7QBpwb5C2HRG2HzhC3RpV3A7DGGNcF6wzjlNVdRwwLkjbizjfrNntdgjGGBMRgnXG8ZyIrBWRx0WkU5C2GRFUlcternS1bsYYU2ZBSRyqOhDPld17gFdFZKWIPBSMbbtt2vIdLNlyAIB7B9tFf8YYE6wzDlR1l6pOAO7Ac03H30p5SYWw73AWAG0b1uTewe1djsYYY9wXrAsATxORR0RkNfAfPD2qmgdj225SVR770jNQ7/S7z3I5GmOMiQzBahx/E/gAGKKqRceaqrCSUw8XLFeNiXYxEmOMiRzBShwDgVOAuiKyX1WPBmm7rko/mg3A05eXOpGhMcacNMpVVSUiMSLyNLAVeAt4F9gmIk+LSIW/6CHtiCdxnNKgpsuRGGNM5ChvG8e/gQSgjar2UtUeeM486gDPlDc4t+1M85w4Na1jEzYZY0y+8iaOEcBtqnoov0BV04E/ABeWc9uuO5Dh6VGVUCPW5UiMMSZylDdxqDN6bdHCXCrBIIf7M7KpERttDePGGOOlvIljjYj8vmihiFwHrC3ntl138EgW8dUqfFONMcYEVXl7Vd0JfCoiNwNL8Jxl9AGqAaPKuW3XHT6aQ604SxzGGOOtXIlDVbcDp4vIeXhGxxXgK1X9LhjBue3Q0RxqxQWrx7IxxlQOwRqr6ntVfVFVJwSaNERkqIisE5FkERlbwnqXi4iKSO/yR+yfw8dyqGmJwxhjCgnaWFVlISLRwERgGNARuFpEOvpYrxaeyaEWhDO+/RlZ1K1uPaqMMcabq4kDzzzlyaq6UVWzgA+BkT7Wexx4GgjbFenZuXnsTDtC87rVwrVLY4ypENxOHM2AbV6PU5yyAiLSA2ihql+GM7CdB4+Sp9CibvVw7tYYYyKe24lDfJQVXP8hIlHAeOBPpW5IZLSILBaRxXv27Cl3YNsOZALQPMHOOIwxxpvbiSMFaOH1uDngPbpuLaAzMEdENgP9gGm+GshV9VVV7a2qvRs0aFDuwLbt9yQOO+MwxpjC3E4ci4B2ItJaRGKBq4Bp+U+qapqq1lfVRFVNBH4BLlbVxaEObNO+DGKihCbxNk6VMcZ4czVxqGoOcBcwE/gN+FhVV4vIYyJysZux/brlIJ2bxRMT7XZuNcaYyOL6RQqqOgOYUaTM57SzqnpuOGLKysljecpBru/XKhy7M8aYCsV+TvuwcvtBjuXk0TuxrtuhGGNMxLHE4cMXy3dSJVo4o019t0MxxpiIY4mjiJQDmby3YAsXdWtKfHUb4NAYY4qyxOHlaHYuj32xhuxc5b7B7d0OxxhjIpLrjeORQlW5edIift64j7HDOtAiwa7fMMYYXyxxOESE2wa04eb+rRncsZHb4RhjTMSyxOFlYIeGbodgjDERz9o4jDHGBMQShzHGmICIqpa+VgUjInuALWV8eX1gbxDDCRaLKzAWV2AsrsBU1rhaqWqpo8RWysRRHiKyWFXDNj2tvyyuwFhcgbG4AnOyx2VVVcYYYwJiicMYY0xALHGc6FW3AyiGxRUYiyswFldgTuq4rI3DGGNMQOyMwxhjTEAscRhjjAmIJQ4vIjJURNaJSLKIjA3zvjeLyEoRWSYii52yBBGZJSJJzn1dp1xEZIIT5woR6RnkWN4QkVQRWeVVFnAsInKDs36SiNwQorgeEZHtznFbJiIXej33oBPXOhG5wKs8aJ+ziLQQkdki8puIrBaRe5xyV49XCXG5eryc7cWJyEIRWe7E9qhT3lpEFjjv/yMRiXXKqzqPk53nE0uLOchxTRKRTV7HrLtTHs6//WgR+VVEvnQeu3qsUFW7edp5ooENQBsgFlgOdAzj/jcD9YuUPQ2MdZbHAk85yxcCXwEC9AMWBDmWs4GewKqyxgIkABud+7rOct0QxPUI8H8+1u3ofIZVgdbOZxsd7M8ZaAL0dJZrAeudfbt6vEqIy9Xj5exLgJrOchVggXMsPgaucspfAf7gLI8BXnGWrwI+KinmEMQ1Cbjcx/rh/Nu/H3gf+NJ57OqxsjOO4/oCyaq6UVWzgA+BkS7HNBJ4y1l+C7jEq/xt9fgFqCMiTYK1U1WdC+wvZywXALNUdb+qHgBmAUNDEFdxRgIfquoxVd0EJOP5jIP6OavqTlVd6iwfAn4DmuHy8SohruKE5Xg58aiqHnYeVnFuCpwHTHHKix6z/GM5BRgkIlJCzMGOqzhh+SxFpDkwHPif81hw+VhZ4jiuGbDN63EKJf+jBZsC34jIEhEZ7ZQ1UtWd4PkiAPKH73Uj1kBjCWeMdzlVBW/kVwm5EZdTLdADzy/ViDleReKCCDheTtXLMiAVzxfrBuCgqub42E9BDM7zaUC9UMRWNC5VzT9m/3CO2XgRqVo0riL7D3ZczwMPAHnO43q4fKwscRwnPsrC2Ve5v6r2BIYBd4rI2SWs63as3oqLJVwxvgycAnQHdgLPuhGXiNQEPgHuVdX0klZ1Oa6IOF6qmquq3YHmeH75nlbCfsIWW9G4RKQz8CDQAeiDp/rpL+GKS0RGAKmqusS7uITth+VYWeI4LgVo4fW4ObAjXDtX1R3OfSrwGZ5/pt35VVDOfaqLsQYaS1hiVNXdzj97HvAax0+/wxaXiFTB8+X8nqp+6hS7frx8xRUJx8ubqh4E5uBpI6gjIvlzBHnvpyAG5/l4PFWWIYvNK66hTrWfquox4E3Ce8z6AxeLyGY81YTn4TkDcfdYlbVxpLLd8ExqtRFPw1F+I2CnMO27BlDLa3k+njrRf1O4gfVpZ3k4hRvlFoYgpkQKN0IHFAueX2ab8DQO1nWWE0IQVxOv5fvw1OMCdKJwY+BGPA29Qf2cnff9NvB8kXJXj1cJcbl6vJx9NQDqOMvVgHnACGAyhRt8xzjLd1K4wffjkmIOQVxNvI7p88CTLv3tn8vxxnF3j1V530xluuHpJbEeT33ruDDut43zoS4HVufvG0/d5HdAknOf4JQLMNGJcyXQO8jxfICnGiMbzy+VW8oSC3Aznka4ZOCmEMX1jrPfFcA0Cn8xjnPiWgcMC8XnDJyF55R/BbDMuV3o9vEqIS5Xj5ezva7Ar04Mq4C/ef0fLHTe/2SgqlMe5zxOdp5vU1rMQY7re+eYrQLe5XjPq7D97TvbPJfjicPVYxXxQ46ISA3gJSALmKOq77kckjHGnNRcaeMQHxdyOeW+LjS6FJiiqrcBF4c9WGOMMYW41Tg+iSL9mkUkGs9p3zA8F6tcLSId8TTi5Hcjyw1jjMYYY3xwJXGo7wu5irvQKAVP8gDrBWaMMa6LKX2VsPF1gcrpwATgPyIyHPiiuBc7F82NBqhRo0avDh06hDBUY4ypfJYsWbJX/ZhzPJISh88LVFQ1A7iptBer6qsishO4qEmTJr0WL14c9ACNMaYyE5Et/qwXSVU/rl6AZ4wxxj/FJg4R+afX8pAwxLIIaOcMFxyL5+KVaYFsQFW/UNXR8fHxZQpgyZYDLNq8nwMZWWV6vTHGnAxKOuPw7vX0VDB3KiIfAD8Dp4pIiojcop4Bue4CZuIZyfNjVV0d4HYvEpFX09LSyhTXc7PW8btXfqbH47MY9dJPLN7s70Csxhhz8ij2AkARWaqeQfcKLVcEvXv31rK0caQcyCQp9TBrdqTz/oKtpB6D0e7oAAAgAElEQVQ6yhd/PIsOjWuHIEpjjIksIrJEVXuXul4JiSMFeA5Po/V9znIBVX3O1+vcJCIXARe1bdv2tqSkpHJt60BGFoOf+4FTGtbk49vPCE6AxhgTwfxNHCVVVb2GZ+awml7L3reIU942Dm91a8QyZmBbFm7aT8/HZ7EyJY3cvMgensUYY8Ih4seqCkQwzzgADh/LYeAzc9hz6BgAXZrFM+2u/ngm1DLGmMolGGccvja6tOwhhV4wzzgAalaNYeFfB/H+racDsHJ7GkPGz2Xm6l1B2b4xxlREgV7HcdL91BYRzmxbn9WPXsA/RnUmJkq4/Z0l3DJpEQczrduuMebkE2jimB6SKIKkvN1xS1KjagzXnt6Kz+/qz/1D2vP9ulSGjJ/L/A17g74vY4yJZJWqjSNfWbvjBmLBxn3c/eGv7E4/xgWdGvHM77pRK65KSPdpjDGhFJI2Dq+NryzL6yqT09vUY9b953D3oHZ8s2Y31/1vAUm7D7kdljHGhFyxgxyKyKXFPQU0Dk045ePVqyos+6sdV4X7h7SnSXwcj32xhguen8tF3Zry2MWdia9uZx/GmMqppAsAs4H38MxbXNTlqhqR13JAeKqqitpx8AiT5m/mf/M2UqNqDMO7NOGibk0585R61n3XGFMhBOPK8SXADaq6ysdz21S1hY+XRQQ3Eke+NTvS+e/cDXy7ZjcZWbn0aFmHy3s1p2uzOrRuUIOaVSNpJHtjjDkuGIljALBFVbf6eK63qkbshBduJo58R7NzmbwkhdfnbWTzvsyC8oa1qtKndQJ9WtWle8u6dG0WT1SUnZEYY9xX7sRRkUVC4sinqmzam8G6XYfYuDeD5NTDzF2/h33O0O3VqkTTtXk8vVrVpWVCdVokVKdF3eo0qRNHlehImi7FGFPZ+Zs4rN4kxESENg1q0qZBzYIyVWV3+jHmJe1h1fY05iXt5dW5G8nxGgsrOkpoEh9Hi7rVaZlQnQ5NatG5WTwdGteybr/GGFdVqjOOYI9VFU45uXnsTDvKtgOZpOw/wrYDmWzdn8m2/Zls2ZdZcIYCcEqDGpzWpDadm8XTqWltOjWNJ6FGrIvRG2MqA6uqipCqqmDIP0NZszONVdvTWZGSxtpd6aQcOFKwTpP4ODo1rU3HpvH0TUygb+sEYmOsqssY47+gVlWJyHmq+n3+ffnDM4EQERrHx9E4Po7zOjQqKD+QkcWanems3pHG6h3prN6RzvdrU8lTqB0Xw/CuTRnSsSH92tSjeqzVShpjgsOvM478GQArykyAle2MIxCZWTn8lLyPL1fs4JvVuzmSnUtsdBS9E+tyTvsGDOzQkHYNa9q1JcaYEwS1qsorcfyqqj2CEmEIncyJw9vR7FwWbz7A3KQ9/LBuD+ucIVGa1anGZT2bcVXfljStU83lKI0xkaLSJA4RaQOMA+JV9XJ/XmOJw7edaUeYs24PX6/axdykPcRECZf3as6dA9vSvG51t8MzxrgspIMcBhDEGyKSKiKripQPFZF1IpIsImNL2oaqblTVW0IZ58miSXw1ru7bkrdu7svcPw/kyj4t+GTJdgY9+wPjZ63nSFau2yEaYyqAUHe7mQQM9S4QkWhgIjAM6AhcLSIdRaSLiHxZ5NYwxPGdtFokVOeJS7ow58/ncn6nxrzwXRKDn/uBGSt3Uhl72hljgsffxHHYuQ9o3HBVnQvsL1LcF0h2ziSygA+Bkaq6UlVHFLml+rsvERktIotFZPGePXsCCfOk1rRONV68ugcfje5HrbgYxry3lGteW8DaXeluh2aMiVB+JQ5VPdv7vpyaAdu8Hqc4ZT6JSD0ReQXoISIPlhDjq8CjwNLYWLsYLlCnt6nHl388i8cv6cxvu9IZPuFHnvtmHdm5eW6HZoyJMG5cIearH2ixdSOquk9V71DVU1T1XyVtWFW/UNXR8fHx5Q7yZBQTHcX1/Vox+0/nMrJ7UyZ8n8wtby22tg9jTCFuJI4UwHtI9ubAjmBsOJRzjp9M6taI5bkruvPkpV34MWkP172+gO0Hj5T+QmPMScGNxLEIaCcirUUkFrgKmOZCHKYUV/VtyYtX92TtznQufGEes9f63eRkjKnESkwcIhItIt+WdeMi8gHwM3CqiKSIyC2qmgPcBcwEfgM+VtXVZd2HN6uqCr7hXZsw454BNKtTjZvfWsTrP26yXlfGnORKvQBQRKYB16tqxNf/VOTRcSNdZlYO9320jJmrd3Ndv5Y8clEnYmy+EGMqlaBdOS4iHwP9gFlARn65qt5d3iBDxa4cD428POWpmWv57w8b6ds6gf9c3YOGtePcDssYEyTBHB13unMzJ7moKOHBYadxaqNajPtsFVe9+guf/OFM6tpcIMacVPwdqyoWaO88XKeq2SGNqoysqip8Fm3ez9Wv/kLX5vG89vve1KtZ1e2QjDHlFLSxqkTkXCAJzzAhLwHrRSQYFwIGnTWOh0+fxAQmXN2D1TvSufTl+RzwmqHQGFO5+dO6+Sxwvqqe41w5fgEwPrRhlY1dxxFeF3Zpwts392XnwaNc/sp8Nuw5XPqLjDEVnj+Jo4qqrst/oKrrgSqhC6ns7Iwj/E5vU4+3b+nLgcxsLnt5Pqu2W9I2prLzJ3EsFpHXReRc5/YasCTUgZmKo1+bekwd058asTFc/dov/Lr1gNshGWNCyJ/E8QdgNXA3cA+wBrgjlEGZiqdlvep8fMcZJNSI5frXF7JkiyUPYyqrUq8cB15X1edU9VJVHaWq41X1WJjiC4i1cbirWZ1qfHBbP+rVjOXq135h5updbodkjAmBEhOHquYCDZzuuBHP2jjc17RONaaO6U+nprW5+4Nf+WG9zY1iTGXjT1XVZuAnEXlYRO7Pv4U4LlOB1a0Ry+s39KFNg5rc+d5S1u0KaP4vY0yE8ydx7AC+dNat5XUzplgJNWL53w29qR4bze/fWMAOG5bdmEqjxCFHnDaOmqr65zDFYyqRZnWq8c4tp3PZy/O55a3FfHx7P2rFRWRPbmNMAPxp4+gZpljKzRrHI8+pjWsx8dqeJO0+xO3vLOFots0maExF509V1TIRmSYi14vIpfm3kEdWBtY4HpnOad+Af/+uKz9v3Mdd7/9Kjs1jbkyF5s/ouAnAPuA8rzIFPg1JRKZSGtWjOYeO5vC3z1fz0NRV/OvSLoj4mn7eGBPpSk0cqnpTOAIxld/vz0hk+4Ej/HfuRlrWq86Yc9u6HZIxpgyKrapyJnDKX36qyHPfhDIoU3mNHdaB4V2b8MzMdcxZZ3OYG1MRldTG0c5reUiR5xqEIBafROQSEXlNRD4XkfPDtV8TGiLCU5d1pUPj2ox5bylrd6W7HZIxJkAlJY6SZngqffYnQETeEJFUEVlVpHyoiKwTkWQRGVvSNlR1qqreBtwIXOnPfk1kq1k1hjdu7EPNqjE2rpUxFVBJiaO6iPQQkV5ANWe5Z/5jP7c/CRjqXeBcGzIRGAZ0BK4WkY4i0kVEvixya+j10oec15lKoHF8HO/eejrVY6O57OX5/HvmWrdDMsb4qaTG8Z3Ac87yLq/l/MelUtW5IpJYpLgvkKyqGwFE5ENgpKr+CxhRdBvi6XrzJPCVqi71Z7+mYmjfqBZTx/Snx+OzmDh7A+0b1WJk92Zuh2WMKUWxiUNVB4Zon82AbV6PU4DTS1j/j8BgIF5E2qrqK75WEpHRwGiAli1bBilUE2p1a8Sy8pHzufWtxdz70TKOZOVyVV/7/IyJZP5cxxFsvjrvF9tmoqoTgAmlbVRVXxWRncBFsbGxvcoRnwmzWnFVmHRTX+54dwljP11JRlYut5zV2u2wjDHF8OfK8WBLAVp4PW6OZyBFcxKrFhvNq7/vxbDOjXn8yzU8+806VP3qg2GMCTM3EscioJ2ItHbm+bgKmBaMDduQIxVb1ZhoXry6ByO7N+XF75MZ895SsnJseBJjIo1ficMZn+o5EXlWREb5u3ER+QD4GThVRFJE5BZVzQHuAmYCvwEfq+rqsgTvY382yGEFFxMdxfgruvOnIe35atUuhk+Yx5It+90OyxjjRUqrDhCRl4C2wAdO0ZXABlW9M8SxlVnv3r118eLFbodhyunzZdu558NlAPxtREdutnYPY0JKRJaoau9S1/MjcawGOquzoohEAStVtVNQIg0iEbkIuKht27a3JSUluR2OCYI9h44x6Nk5pB/N4cIujXnuiu7EVYl2OyxjKiV/E4c/VVXrAO/+kS2AFWUNLJSsjaPyaVCrKvMfHMSlPZoxY+UuLnrxR5uK1hiX+ZM46gG/icgcEZkDrAEaOHN0BKVRO1isjaNyqlk1hueu7M5fL+xAUuphLnh+Lsu3HXQ7LGNOWv5UVZ1T0vOq+kNQIwoCa+OovP43byNPTP8NgH+M6sy1p7dyOSJjKo+gtXFURJY4Kre1u9IZ+vw8AFomVOeb+862dg9jgiBobRwi0k9EFonIYRHJEpFcEYnIsbCtqurk0KFxbdY+PpRereqydX8mHR7+mpQDmW6HZcxJw582jv8AVwNJeEbFvdUpizjWOH7yiKsSzZQ7zmD02W0AOOup2Xy0aKtdbW5MGPh1AaCqJgPRqpqrqm8C54Y0KmP8ICL89cLTmHH3AAD+8slKrnt9Acdycl2OzJjKzZ/EkekMDbJMRJ4WkfuAGiGOq0ysqurk1LGpp+qqfaOa/JS8j1Mf+preT3zLvsPH7AzEmBDwJ3FcD0TjGSYkA891HJeFMqiysqqqk1dclWhm3ns2Dww9FYC9h4/R64lv+fu0oIxmY4zxYr2qTKWzfvchzh8/t+Bxi4RqvH9rP1okVHcxKmMiX7m744rISkqeJ6Nr2cMLLUscBuDBT1fwwcLjc4bVjoth+d/PxzOppDGmqGAkjhKvrFLVLWWMLeQscZh8x3Jyufzln1m5/Xi718BTG/DmTX1djMqYyBSUCwBFJBqYqaqDgxlcqNggh6Y42bl5tBv31Qnlb97Yh4EdGroQkTGRJ5ij404DrlfVCtNVyc44THHSjmRz1lPfc+hoTqHyTk1r8/md/YmJdmNuM2MiQzATx8dAP2AWnl5VAKjq3eUNMlQscZjSbN2Xydn/nu3zuWd/143E+jXo1apumKMyxl3BTBw3+CpX1bfKGFvIWeIw/lJVvlmzm9vfWeLz+Q6Na/HPS7tw6UvzAUj+x7ATzkqycvKIjQnOmYqqkpunEXvmo6os2LSf01snWCcDP63dlU6UCO0b1SrzNvK/p0N9zIM6yKGIVANaquq6YAQXapY4TFnk5Sk3TVrED+v3lLjeJd2bcv+QU2lZrzozVu5kzHtLeeeWvvRJTOCvn65kaOfGnN+pccH6yamHiImKIrF+6dfNPjR1Je/+spUre7fgqcsjr+PilCUp/N/k5Yy/shujejR3O5wKIXHsdAA2Pzm8zNvo8shMGteOY9b9JQ5WXm7BHOTwImAZ8LXzuHukzcNhTDBERQlv3dyXzU8OZ90TQ3n0Yt+TXE5dtoOz/z2bxLHTGfPeUgCuf30hHR7+mk9/3c7od5bw+o+bUFVempPM4Ofmcu4zc8jLO/FHWm6eFgyRcvcHv/LuL1sB+GjxthPWDaWUA5lc/eovpB3JLnad+Rv28n+TlwOwZV/oBpX8fNl22j/0FUezI3fomNw8ZVkZ5oRJTT/Kewv865C659AxXvg2ibw85dDRHJJSDwe8v1Dx53z4EaAvcBBAVZcBYZv8WUROE5FXRGSKiPwhXPs1J7eqMdHccGYim58czvonhrH5yeF0alrb79c//uUaWj84g6e/Pn6S/v3aVJZuPVDw+Gh2Lo9MW82pD31N2pFspi3f4XNbObl5pGUW/kLftDeDB6YsJyc3j3/PXMu3a3YDni+mNTv8G7z6SFYuHy/ehqry/LdJ/LxxH1+v2llonc9+TWGrkySueW1BQbmPHBg0T0z/jaycPPZnZBUqf/eXLQx9fi47044U+1pVZfLibWQcyyl2naJSDx0lcez0Yr/Qk1MP86vX5wbwwrfruWTiT6xMKbnPUNEanb7//I5xn61i2/7SE+9fPlnB+G/X88umfaWuG27+JI4cHz2q/PqzEZE3RCRVRFYVKR8qIutEJFlExpa0DVX9TVXvAK4ASj2FMibY8tsvpt89gM1PDmfzk8OZP/Y87h/SPqDt3Pr2Yi59aT4PTFnOW/M30+Hhr3nnF8+XVbdHvzlh/eXbDpKdm8eIF3+k22PfMPaT4zM2D3xmDh8vTmHl9jQmzt7ArW8vJjX9KH3/+R0XTvDMVXLn+0tJHDudtbvSyctT/vTxcqb+ur1gG6f97WsemLKCqcu2FySt9xdsJS9PUVUOZmZx30fLefQLH8O2BDjixO70o0xZklLoizQtM9vnWGJ7Dh0D4Eh2LrvSjifCh6auYu2uQ5zxr+/5ZaPvL9P3Fmzlz1NW0OnvM8n1ym5Juw8VO+Xw6u2e7Y/7bBXPzFzHxf/5sdDzg5/7gVFOG1e+RZs9iWR/ZuHk5i03T/kxea/P5+YlHS9/5YcNzFmXesI636/1lHkn7OJ8s3oX/5u3sdT1giXGj3VWicg1QLSItAPuBuaX8pp8k/AMwf52foFzbchEYAiQAixyqr6igX8Vef3NqpoqIhcDY4nQ4dzNyadpnWrcPagddw9qV1CWnZtHlegotu3PZMDTvntsAXy8OAXPn37JRk78qdDjDxdt457B7Xh/wdaCsmSv6ou+//yuYHl/RhbTV3jOHvInvQL4ZGkKgzs2Itar8f2+j5YXLC9PSaPnE7NoVa8G+c2w36098UvNnzOOLfsymDg7mUcv7szpTmyt6lWnT2JCQa+2ewe3497B7Qvey+Dnjk8oevhoDoOe9Tz+8o9nFdr29BU76dem3gn7/HnD8YTy5YodjOzeDIAhzhA0vtoZ0o8eP5v7z+xkoPR2iZ+dxPX6j5to27Ama3akM6Rjo0LrXDLxp0IXnnqfKT3+5RquOb0la3ak8+RXa0vcV2myc/MY7XTuuHVAmzJtI1D+JI4/AuOAY8D7wEzgCX82rqpzRSSxSHFfIFlVNwKIyIfASFX9FzCimO1MA6aJyHQnhhOIyGhgNEDLli39Cc+YoKrifBm3SKhe6EtAVWn94Ay/tvHnC07lvz9sIP2o76qWM/71feH1p6zwuV7Px2cVu4/Of59J/ZqxxT5/MDObg5mF6+8PF6n6ySvljCPjWA6v/LCRjxenFOrW/NGibdzwxkL6t60PwPPfJhUkDu+kAbAi5XgMRdteDh/LIXHsdJ66rAtX9jn+/z595fGqtvQS2mu8RUcV31PJu4rQV++5uev3cOlLP7E7/Rib/nVhoV5P3kkDYNX241WIg07zXHS6zccEZEeycnn481UnlBfnL8X8DYSSP4njVFUdhyd5BEMzwLvlLwU4vbiVReRc4FKgKlDsf5+qvioiO4GLYmNjewUnVGPKT0R8/ppMy8wmJlqoViWa299dQuem8dw5sC19Wyfwu1d+DmlMew8XX8XiS+e/zyz0OP8CyqPZuT6n7e3ktf5fPllZsDxliedMa5bTJlOShz8/XkV27f8KV9fsTj8KwMTZGwolDm+ZWbmM/WQFbRvWLChLHDv9hC7VG1IzfL0cgD7//LZgeenWA1SJFi57ufBnszvdU7W2bf8RWtYrfiDN294+3tPz1Ea1SDmQWagb+IY9h3nqq7Us2XKAfRklfz55eUquKlWio/jUq/oxXPxJHM+JSBNgMvChqpZ3nGpf6b2kwRTnAHP82bCqfgF80bt379vKFJkxYRRfvUrB8mu/P9581ycxgc1PDmd/RhaHj+Ywe10qUVHCw1OL/xV618C2TF22nZQDhRuOa8fFFHv2Uh4Na1Xlvo+W8dmv25n3wEBaJFRn0LNz2LAngyl3nBHQtvo/+T3bDxbf4O1L/llCdm5esev8y6kCKuqDhVuZvnIn7RrWIjpKmDR/c7HbyMo5vv2qMVF8sHBrsevmX1A65txTeGnOhpLC59lZ63l21vpCZWPeXcq63b7bYYr6w3tLWLBpPweLdJpQ1bBcX1Nq4lDVgSLSGE/j9KsiUhv4SFX9qq7yIQXPnB75mgO+u5MEyGusqmBszhhXJdSIJaFGLDecmQjAyO5N2XHwCEeychn10nyGdmpMn9YJVIkWfn9GIv93gWcukpUpaVzkNPCueOQC3v1lCw9NXcVTl3Up9Osf4IJOjZi5uvRf//keHtGRx79cQ6PacQVffAOens0VvZuzYY/nl/vlAZ4tBZo0gIIvTO8v9vxG9dLkn8n8snF/QPvMzlW/jlVpSaM4pSWNpvFxBcvFxXEsJ8/nGWCw+XPGgaruAiaIyGzgAeBv+NnO4cMioJ2ItAa2A1cB15RxW8acNGrHVaF2Y89ZSkkNqV2axxd6/rp+rbiun2ew60t6NOPUh76mR8s6XNS1KTef5elZn+P8cr/mtQX8ZVgHLnvZ0/9l5r1nczQ7t6ChflSPZjz+5Roe+KRwvbqnwT988tsPsrzOOPr849viVg+KzKycEq9z8UeHxrVYW0zvLl8+HXNmwagFqYeOsXlvBnWrF98+NfXX7VzVN/RtvP4MOXIacCXwO2Av8CHwiaqe2NXixNd+gGd+8vrAbuDvqvq6iFwIPI+nJ9UbqvqP8ryJouzKcWPKJy9Pyco9/ut1RcpBWtStTt0asQU9jgL10PDTeGL6bwWPY6OjCn3xA3xz39nUjqtCv399V/Tl3Ny/NTeemXjCGGObnxxOTm4ebZ3Rjzs2qc2anf5dy+LLonGD2Z1+lEsm/kROkC9Yubhb02Kv1ykulk+XphRb7eZLea5QD9qV48CbwAFgiKqeo6ov+5M0AFT1alVtoqpVVLW5qr7ulM9Q1faqekowk4bNOW5McERFSaEqj67N61C3hueX7vu3Hu/L8tPY8/jjeb6rhr/70zks//v5gKet5dYBbdjwzwsLnr/j3FMA+OuFHQBPl9v2jWrROD6OPzvVbq9c14uF4wbx+MhO/O2ijjSoVfWE/bz7y5aCpAHw/m3F9rUp1ps39ilYrl8zls7N4pnyhzMD3o4vl/c6PjRLrbgYPhrdz+/X1qsRy+3nnFLqet4dw3yNUBBs/pxxVANOwdOAvUFVj4Y8qnKyMw5j3LF+9yEmL97GsC5N6Nmy5NGF8/KUrfsz/RrDy1tJZzwvXduTC7s04cIX5hETLUSJFAwNsvzv5xdcaNmgVlUWjRtcsK2Z955N4/g4jmXn0rD28baEfYePkZ2rJ5wBfTS6H6c1rc1V//2FNTvTaVO/Bhv3ntg7a/OTwwt1xx7UoSGv39iHb9fs5lanl9XlvZoX9DbLN+mmPpx76vF5Yh77Yg1v/LTJ53uuHhvNy9f14oY3FgLw7f1n07Zh2QZU9PeMo9g2DhGJAf4J3ARsxXN20lxE3gTGqWr5KvuMMZVO+0a1GDe8o1/rRkVJwEkDYOI1Pflo8Tbm+hiM8sIuTQBP20BMlJBxLJduj3mSRe24GNY9MZTNezM5pUHh/Z7a2PmirValUHm9mlXJzCrcK23RuMEFZz6fjjmTzKxcEmrE8sP6Pdz61iKycz0/xl+5rifg6Y49+uw2vDp3Iz2da1oGd2zE5ieHF3Rnzk8cH99+BvM37C2UNADGDutwQuI4rUlt/jSkPYOLXHgYjl5VJU0dOx6oBdynqoecstrAM8ARVb0n5NEFyGYANObksuPgEf49cx0HM7N4eERH2jSoecI6c9al0s2rqs3b+t2H2Hc4izNOOfEqdG9rdqQTHSXHE0wxlmw5wGUvz+e+we25Z3DhUQXmb9jH2e3q+/xif2lOMk3jq3FJj2Ylbv/88T+wfvdh7h7U7oQhb579Zh0vfp98woWIgQjGnONJQHstsoIzZMhaVW3n84URwKqqjDEmcMFoHNeiScMpzMXPQQ7DzRrHjTEm9EpKHGtE5PdFC0XkOsD/vmFhpKpfqOro+Ph4t0MxxphKq6QLAO8EPhWRm4EleM4y+gDVgFFhiM0YY0wE8qc77nlAJzxjTK1W1ROvzIkwIrIH8G+arRPVx3OhY6SxuAJjcQXG4gpMZY2rlao2KG0lv+YcP5mIyGJ/GofCzeIKjMUVGIsrMCd7XP5cOW6MMcYUsMRhjDEmIJY4TvSq2wEUw+IKjMUVGIsrMCd1XNbGYYwxJiB2xmGMMSYgljiMMcYExBKHFxEZKiLrRCRZRMaGed+bRWSliCwTkcVOWYKIzBKRJOe+rlMuIjLBiXOFiPQMcixviEiqiKzyKgs4FhG5wVk/SURuCFFcj4jIdue4LXMmCct/7kEnrnUicoFXedA+ZxFpISKzReQ3EVktIvc45a4erxLicvV4OduLE5GFIrLcie1Rp7y1iCxw3v9HIhLrlFd1Hic7zyeWFnOQ45okIpu8jll3pzycf/vRIvKriHzpPHb1WKGqdvO080QDG4A2QCywHOgYxv1vBuoXKXsaGOssjwWecpYvBL7Cc1FmP2BBkGM5G+gJrCprLEACsNG5r+ss1w1BXI8A/+dj3Y7OZ1gVaO18ttHB/pyBJkBPZ7kWsN7Zt6vHq4S4XD1ezr4EqOksVwEWOMfiY+Aqp/wV4A/O8hjgFWf5KuCjkmIOQVyTgMt9rB/Ov/37gfeBL53Hrh4rO+M4ri+QrKobVTULzxS5I12OaSTwlrP8FnCJV/nb6vELUEdEmgRrp6o6F9hfzlguAGap6n5VPQDMAoaGIK7ijAQ+VNVjqroJSMbzGQf1c1bVnaq61Fk+BPwGNMPl41VCXMUJy/Fy4lFVPew8rOLcFDgPmOKUFz1m+cdyCjBIRKSEmIMdV3HC8lmKSHNgOPA/57Hg8rGyxHFcM2Cb1+MUSv5HCzYFvhGRJSIy2ilrpKo7wfNFAOTP7uJGrIHGEs4Y73KqCt7IrxJyIy6nWqAHnl+qEXO8isQFEXC8nKqXZUAqni/WDcBBVc2fNcl7PwUxOM+nAfVCEVvRuFQ1/5j9w4phkHsAAAUFSURBVDlm40Ukf/7acB2z54EHgPwJ2uvh8rGyxHGcr5lPwtlXub+q9gSGAXeKyNklrOt2rN6KiyVcMb6MZ2rj7sBO4Fk34hKRmsAnwL2qml7Sqi7HFRHHS1VzVbU70BzPL9/TSthP2GIrGpeIdAYeBDrgGeQ1AfhLuOISkRFAqqou8S4uYfthOVaWOI5LAVp4PW4O7AjXzlV1h3OfCnyG559pd34VlHOf6mKsgcYSlhhVdbfzz54HvMbx0++wxSUiVfB8Ob+nqp86xa4fL19xRcLx8qaqB4E5eNoI6ohnyuqi+ymIwXk+Hk+VZchi84prqFPtp6p6DHiT8B6z/sDFIrIZTzXheXjOQNw9VmVtHKlsNzxDzG/E03CU3wjYKUz7rgHU8lqej6dO9N8UbmB92lkeTuFGuYUhiCmRwo3QAcWC55fZJjyNg3Wd5YQQxNXEa/k+PPW44BnR2bsxcCOeht6gfs7O+34beL5IuavHq4S4XD1ezr4aAHWc5WrAPGAEMJnCDb5jnOU7Kdzg+3FJMYcgriZex/R54EmX/vbP5XjjuLvHqrxvpjLd8PSSWI+nvnVcGPfbxvlQlwOr8/eNp27yOyDJuU9wygWY6MS5Eugd5Hg+wFONkY3nl8otZYkFuBlPI1wy8P/t3cGLVWUYx/HvD4VSgyTSbSJMthBxF4qjs4gkkTYSQougTeYiQSgXCW0MFP0LbFsaolAJgm5CxcLGKOcySUFLEVuki1SKkqfF8xycO17n9srcZuHvA8PMPfee97xzZu557nvPPb/37RH169Pabg84Tf+BcX/16xfgtVH8nYFN5JC/B1ytr20Lvb/m6NeC7q9qbx3wY/VhGvhoxvNgsn7/k8BTtfzpuv1r3b96WJ/nuV9f1z6bBj7jwSev/rf//WpzggeFY0H3lSNHzMysic9xmJlZExcOMzNr4sJhZmZNXDjMzKyJC4eZmTVx4TADJN2p76skvTnPbX846/a389z+mkpw1Xy3bTaIC4dZv1VAU+GQtGjIQ/oKR0RsbOzTMOPkxWrryOuAzEbKhcOs3yFgvOZd2Fuhd0ckXamQu10AkiaU810cJy/+QtKXFVL5UxdUKekQsKTaO1bLutGNqu1p5VwsO2e0fV7SKUk/SzpWCad9JI1XIN9h4H3gDLBVNZ+L2aj4AkAz8mAeEc9ImiDnq9hey98BVkbEx5WK+g3wBvACeaBeGxlTjaTnIuKWpCXAFWBLRPzetT1gWzuAd8l4medrnZeBNcBXZEzEjdrmBxFx6RF9vwxsIHOUjkSERx02Uh5xmM3tVeCtemX/HRklMlb3TXZFo+yRNAVcJgPlxpjbJuDzyNDB34ALZAJr1/b1yDDCq+RbaA+RtBT4M/IV4BgZJ2E2UouHP8TsiSbgvYg417cwRyZ3Z91+BdgQEfcknSdzg4a1/Sh/zfj5PgOeq5JOk3HfyyX1yOLyvaSDEXFiyLbNHptHHGb9/iCnWu2cA3ZXRDmSXpS0bMB6zwK3q2i8RKaldv7u1p/lIrCzzqOsIKfGnfyvHY2I18lo9N3AHjIVdb2Lho2aC4dZvx7wj6QpSXvJ6TqvAT9ImgaOMnikfhZYXK/8D5BvV3U+AXrdyfEZvqjtTZEJrPsi4mZjfzcDl8hPVl1oXNfssfjkuJmZNfGIw8zMmrhwmJlZExcOMzNr4sJhZmZNXDjMzKyJC4eZmTVx4TAzsyb/AuQSDZP64cb+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "J2 = 0.4\n",
    "\n",
    "plt.subplot(2, 1, 1)\n",
    "plt.title(r'$J_1 J_2$ model, $J_2=' + str(J2) + '$')\n",
    "plt.ylabel('Overlap = F')\n",
    "plt.xlabel('Iteration #')\n",
    "\n",
    "plt.plot(iters, overlap)\n",
    "plt.axhline(y=1, xmin=0, xmax=iters[-1], linewidth=2, color='k',label='max accuracy = 1')\n",
    "\n",
    "plt.legend(frameon=False)\n",
    "\n",
    "plt.subplot(2, 1, 2)\n",
    "plt.ylabel('Overlap Error = 1-F')\n",
    "plt.xlabel('Iteration #')\n",
    "plt.semilogy(iters, 1.-overlap)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result suggests that indeed we could have a good approximate to the state given by supervised learning. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
